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Abstract 


In  this  project,  the  main  objective  was  to  develop  multiresolution  wavelet  algo¬ 
rithm  to  study  the  flame  acceleration  and  to  understand  the  mechanism  of  the  tran¬ 
sition  of  deflagration  to  detonation.  In  the  first  half  year  of  the  project,  The  author 
completed  the  study  of  two  dimensional  detonation  waves  based  on  hybrid  high  order 
methods  (a  combination  of  high  order  Essentially  Non-osciUatory  (ENO)  methods  and 
spectral  methods  and  Shock  Tracking  methods).  The  result  on  the  detonation  waves 
was  published  in  the  AIAA  Journal.  In  the  remaining  time  of  the  two  years,  first, 
the  author  completed  the  theoretical  and  algorithmic  studies  of  the  adaptive  wavelet 
method,  which  could  handle  nonperodic  boundary  conditions  and  nonlinear  time  de¬ 
pendent  PDE’s  (The  result  was  published  in  the  SIAM  Journal  of  Numerical  Analysis). 
Secondly,  the  author  implemented  the  adaptive  wavelet  methods  for  the  solution  of  one¬ 
dimensional  flame  propagation;  thirdly,  the  author  developed  a  Fortran  code  WL2D 
(more  than  13,000  lines)  for  the  two  dimensional  multi-scale  wavelet  algorithms  with 
an  efficient  data  structure  and  implemented  a  second  order  implicit  factorized  scheme 
for  the  adaptive  wavelet  methods. 
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1  Hybrid  High  Order  Methods  for  Detonation  Waves 

Partially  supported  by  this  grant,  we  completed  the  development  of  high  order  hybrid  nu¬ 
merical  simulation  of  two  dimensional  detonation  waves.  The  major  finding  of  this  work  was 
that  the  cellular  structure  of  detonation  waves  depended  very  sensitively  on  the  numerical 
dissipation  of  the  algorithms  representing  detonation  fronts.  Further  studies  of  the  work  in 
three  dimensional  cases  were  needed  to  understand  the  three  dimensional  effects  of  detona¬ 
tion  waves.  One  paper  summarizing  the  results  of  this  study  was  published  in  AIAA  Journal, 
Vol.  33,  Number  3,  pp  1248-1255. 

2  Parallel  Multi-scale  Wavelet  Algorithms 

In  an  attempt  to  design  multiscale  methods  for  the  study  of  deflagration  to  detonation  tran¬ 
sition  (DDT)  problem,  we  constructed  a  wavelet  collocation  multi-resolution  algorithm  for 
the  initial  value  boundary  problem  of  nonlinear  PDF’s.  The  key  component  in  this  colloca¬ 
tion  method  was  a  so-called  “Discrete  Wavelet  Transform”  (DWT)  which  maped  a  solution 
between  the  physical  space  and  the  wavelet  coefficient  space.  The  DWT  transformation  only 
took  0{NlogN)  operations  where  N  was  the  total  number  of  unknowns.  Therefore,  the  non¬ 
linear  term  in  the  PDF  could  be  easily  treated  in  the  physical  space,  and  the  derivatives  of 
those  nonlinear  terms  then  computed  in  the  wavelet  space.  The  wavelet  collocation  methods 
had  the  following  advantages:  (a)  the  capability  to  handle  arbitrary  non-periodic  boundary 
conditions;  (b)  the  capability  to  treat  general  nonlinearity  through  the  collocations  of  the 
PDF’s;  (c)  flexibility  of  adaptive  meshing  in  regions  where  high  gradients  of  solution  occur 
as  in  shock  waves  and  turbulent  premixed  flames;  (c)  Capability  of  parallel  and  multi-grid 
implementations  . 

The  following  paragraph  summarizes  the  key  technical  details  of  the  algorithms.  One 
paper  containing  some  of  the  details  of  the  results  of  this  study  was  published  in  the  Journal 
of  SIAM  Numerical  Analysis,  Volume  33,  Number  3,  pp.  937-970,  June  1996. 

2.1  Discrete  Wavelet  Transform  (DWT) 

The  Discrete  Wavelet  Transformation  (DWT)  maps  between  a  function  at  its  sample  points 
to  the  coefficients  of  its  wavelet  interpolation  expansion.  This  transform  in  both  direction 
will  only  take  0{NlogN),  N  =  +  L  —  1  operations  for  the  77^(7) -wavelet  basis. 

Let  f{x)  in  Hq{I),  and  we  intend  to  construct  a  wavelet  interpolation 

Vjfix)  €  Vo  ®  Wo  ©  •  •  •  ©  Wj  for  J  >  0 

First  introduce  two  index  spaces,  for  Wj,j  >  0  we  define 

X:,  =  {-l,0,---,n,-2}  (1) 

where  nj  =  DimWj  =  2^  L  and  for  Vo  we  define 

/C_i  =  {-l,0,---,L-3}  (2) 


3 


Second,  define  the  collocation  points  in  domain  /  =  [0,  L]  associated  with  wavelet  space 


VEi,i  >  0  by 

(i)  ^  +  1-5 

^k  -  ’ 

k  G  }Cj 

(3) 

and  for  Vb  the  collocation  points  are 

=  k  +  2, 

k  G  /C_i 

(4) 

The  interpolation  operator  ly^f  in  Vb  interpolates  f{x)  on  collocation  points 
while  the  interpolation  operator  1^,^/  in  Wj  interpolates  any  function  f{x)  in  Hq{I)  on 

collocation  points  The  construction  of  the  interpolant  and  Iv^f  involves  the 

inversion  of  a  tridiagonal  matrix  with  size  dimWj  =  nj  or  dimVo  =  L  —  1,  respectively. 

Now  let  us  assume  that  the  values  of  a  function  f{x)  6  Hq{I)  are  given  on  all  the 
collocation  points  G  <  j  <  J,  then  the  wavelet  interpolation  Vjf{x)  G 

Vb  0  ITo  0  •  •  •  0  Wj  for  J  >  0  is  defined  as 

'Pjfix)  =  f-i{x)  +  fj{x)  (5) 

j=0 

such  that 

=  fi^k^)^  for  k  G  JCj, -1  <j  <J 

where 

f-i{x)  =  lvof{x)  =  Y  f-i,k4>o,k{x)  G  Vb 

fcex:_i 

and  for  i  >  0, 

M^)  =  !»,(/'''  -  (Pi-./)''')  =  E  e  Wj.  (6) 

keKlj 

Let  us  denote  f  =  the  values  of  f{x)  on  all  interpolation  points  on 

all  levels,  =  {f{x^k^)}keiCj,j  >  —1  and  f  =  the  wavelet  coefficients 

in  the  expansion  (5),  =  {fj,k}keiCj,j  >  -!• 

A  fast  discrete  wavelet  transform  (DWT)  was  proposed  for  an  efficient  transformation 
between  point  values  f  and  f  with  operation  counts  of  order  O(A^lnA^)  in  both  directions  . 
The  efficiency  of  the  transformation  is  based  on  the  fact  that  (j).i^k{x),'ipj^k{x),k  G  ICj  form 
a  hierarchical  nodal  basis  on  all  levels  of  collocation  points  —l‘^j<J,k^}Cj. 

2.2  Wavelet  Collocation  Methods  for  PDE’s 

We  consider  a  collocation  method  based  on  the  DWT  transform  for  time  dependent  PDE’s. 
Let  u  =  u{x,t)  be  the  solution  of  the  following  initial  boundary  value  problem 

{Ut  0  —  Uxx  0  G  [0,  L j ,  f  ^  0 

ti(0,f)  =go{t) 

Bu{L,t)  =g,{t) 

Bu{x,0)  =  f{x) 
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where  B  is  the  boundary  condition  operator  which  could  be  either  Dirichlet  boundary  con¬ 
ditions  or  Neuman  type  or  Robin  type  boundary  conditions. 

The  numerical  solution  uj[x,t)  will  be  represented  by  a  unique  decomposition  in  the 
cubic  spline  wavelet  decomposition  of  H^{I)  -  14  0  ITo  0  •  •  •  0  Wj,  J  >  0,  namely. 

uj{x)  =  lb,ju{x)  +  ii_i(ar)  +  uo{x)  + - h  uj{x)  (8) 

where  cubic  spline  consists  of  the  nonhomogenuity  of  u{x,t)  on  both  boundaries, 

and  the  coefficients  Uj^k{t)  are  all  functions  of  t.  Using  the  DWT  transform,  we  can  also 
identify  the  numerical  solution  uji^x^t')  by  its  point  values  on  all  collocation  points,  we  put 
all  these  values  in  vector  u  =  u(t),  i.e. 

u  =  u(t)  =  •  •  • , 

where  =  {u{x^,^\t)},k  e  fCj  for  j  >  -1. 

To  solve  for  the  unknown  solution  vector  u(t),  we  collocate  the  PDE  (7)  on  all  collocation 
points,  then  we  have  the  following  semi-discretized  wavelet  collocation  method. 

Semi-Discretized  Wavelet  Collocation  Methods 


'  ujt  +  fx{uj) 

Buj{0,t) 

Buj{L,t) 


=  ^o(^) 
=  9i{t) 


(9) 


where  k^fCj  —  l<j<J. 

Computation  of  fx{Xk^)  =  fx{uj{xk^) 

Step  1  Given  u  =  u(°),  •  •  ■ ,  compute  =  {/(u?)},  k  €  KjJ  >  -1  and 

define 

f  =(f(-l)j(0),...,f(J))T; 


Step  2  Compute  the  wavelet  interpolation  expansion  using  DWT  transform  for  f  as  in  (5) 

Step  3  Differentiate  the  interpolation  expansion  and  evaluate  at  all  collocation  points  x]^\ 

which  is  taken  as  fx{{x[^^))  =  fx{uj{x[^^)). 

The  total  cost  of  computing  the  derivatives  will  be  (5J  +  12)77  <  5NlogN. 


2.3  Adaptive  Meshing 

In  equations  (8),  uj(x)  is  expressed  using  the  full  set  of  collocation  points  As  most 

of  the  wavelet  expansion  coefficients  Uj  for  large  j  can  be  ignored  within  a  given  tolerance 
c.  So  we  can  dynamically  adjust  the  number  and  locations  of  the  collocation  points  used  in 
the  wavelet  expansions,  reducing  significantly  the  cost  of  the  scheme  while  providing  enough 
resolution  in  the  regions  where  the  solution  varies  significantly. 
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Let  e  >  0  be  a  prescribed  tolerance  and  i  >  0,  ^  =  £(c)  =  mm(^,  —  log  e/ log  a),  a  = 
13.98. 

Step  1.  First  locate  the  range  for  the  index  k, 

-  =  m(j,e)  (10) 

such  that 

\uj,k\<^,  k'i  <  k  <  I'i,  i  =  1,- ,m. 

Step  2.  Ignore  Uj^k  in  (8)  for  fc,  <  k  <  li,i  =  !,•••  ,m,ki  =  k'-  +  £  +  3Ji 
namely  we  redefine  Uj{x)  as 

k€ic,\icr 

where  K,'-  =  Ui<i<m[^i,  ^i]- 

Step  3.  The  new  collocation  points  and  unknowns  will  be 

e  /C_i  if  j  =  -1;A;  e  if  j  >  0. 

2.4  Data  Structure  for  Two  Dimensional  Adaptive  Wavelet  Algo¬ 
rithms 

We  have  developed  an  efficient  data  structure  for  two  and  three  dimensional  wavelet  ap¬ 
proximations.  Being  an  adaptive  scheme  in  nature,  an  efficient  data  structure  for  handling 
numerical  solutions  in  a  wavelet  framework  is  the  first  step  in  making  a  theoretical  poten¬ 
tial  into  real  application.  We  have  used  sparse  matrix  data  structure  (compressed  row  and 
compressed  column  vector  technique)  in  treating  the  data  structure  on  each  wavelet  space 
Wf  X  w]  in  the  case  of  two  dimensional  approximations.  Such  approach  have  the  advantage 
of  only  storing  the  mesh  points  used  in  the  adaptive  wavelet  approximations  and  easy  access 
to  numerical  data  on  each  constant  x  and  y  lines  for  the  numerical  differentiations. 

Data  Structure  -  1-D 

One  Dimensional  Case: 


(11) 

=  l'i-£-3, 


u  =  {u^ 

GFo0Wo0Wx0---Wj 


where  nj  -  number  of  mesh  points  on  level  Wj 

[1]  pointer{j)  -  pointer  of  first  element  of 

[2]  npts{j)  -  number  of  elements  in  —  Uj 

[3]  index{l  :  rij)  -  collocation  point  location  indices 
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Data  Structure  -  2-D 


xy-^  , 

Mesh  on  level  {jx,jy)  -  1 


'y  =  .. . 

jy  jx 

JX  ^  jy 


X 

jy 


•  Compressed  Row  Form  : 


i-th  row: 

lenrow{i,jx,jy)  -  length  of  i-th  row 
ipr{i,jx,jy)  -  pointer  of  1st  mesh  on  i-th  row 
icn{:,jx,jy)  -  column  indices  of  i-th  row 
•  Compressed  column  Form  : 
i-th  column: 

lencol{i,jx,jy)  -  length  of  i-th  row 
ipc{i,jx,jy)  -  pointer  of  1st  mesh  on  i-th  row 
irn{:,jx,jy)  -  column  indices  of  i-th  row 


2.5  Fast  Time  Integration  -  Factorized  ADI  Wavelet  Approach 

On  the  issue  of  time  integration,  we  have  successfully  implemented  a  second  order  implicit 
factorized  scheme  of  Beam  and  Warming  type  for  the  adaptive  wavelet  methods.  Unlike  many 
other  adaptive  methods  (finite  element  and  finite  difference),  the  adaptive  wavelet  methods 
actually  can  be  easily  implemented  using  an  ADI  approach.  So,  solution  of  the  algebraic 
systems  from  the  implicit  discretization  of  2-dimensional  diffusion  operators  is  replaced  by 
that  of  only  1-dimensional  operators,  thus  speeding  up  the  time  integrations  tremendously. 
Example:  Beam  Warming  schemes  for  Euler  Equations. 


^  dFju)  dG{u) 
dt  dx  dy 


0 
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3  Computational  Results  of  Multi-scale  Wavelet  Algo¬ 
rithms  for  2-D  flames 

We  have  developed  a  two  dimensional  code  WL2D  based  on  the  adaptive  wavelet  approx¬ 
imations.  The  following  eight  pages  summarize  the  simulation  results  of  two  dimensional 
cellular  and  planner  flame  propagations. 

[1]  Two  Dimensional  Perturbed  Cellular  Flame 


Qt 

Ct 


A0  —  K - h  D 

a 

- K - It 

L  r 


where  A  is  the  Laplace  Operator  Initial  Condition: 


0  =  (^r+o(F) 

c  =  (l-0)  +  O(i) 


+  perturbation  of  flame  fronts 


I 


[2]  Two  dimensional  Perturbed  planar  flames 


Qt 


Ct 


“(■  n 


where  A  is  the  Laplace  Operator  Initial  Condition: 


0  =  exp  a:  if  a:  <  0 
=  1  if  a;  >  0 
C  =  1-0 
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3,4,5  Reactant  Mass  Fraction 
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6  Appendix:  2-D  wavelet  Computer  Code  WL2D 
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